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Following a cell-method due to van Kampen for the calculation of a coarse-grained 
free energy functional for the van der Waals gas, we compute a corresponding entropy 
functional from microscopic principles. This entropy functional is one of the building 
blocks of the recently developed GENERIC framework. This framework allows to obtain 
in a thermodynamically consistent way the continuum hydrodynamic equations for a 
fluid able to display liquid- vapor coexistence. Surface tension appears naturally and the 
resulting model describes interfaces as diffuse regions, much in the same spirit as the 
gradient theory for equilibrium situations. We suggest that using interfacial forces in the 
integral form obtained in the microscopic derivation instead of third order derivatives of 
the density field might represent an advantage from a computational point of view. 



I. INTRODUCTION 

Boiling of water in a pot, the formation of a cloud, or the rise of bubbles in a pint of beer are fascinating 
everyday experiences that involve liquid-vapor coexistence in non-trivial flow situations. Knowledge of 
the dynamics of this phenomena is crucial in less pleasant situations like the prediction and control of 
nuclear accidents in refrigerated nuclear reactors. The complexity of the problem requires the aid of 
computer simulations in order to extract practical information. Engineering conventional computational 
fluid mechanics approaches have relied on effective hydrodynamic equations in which the presence of 
bubbles or drops is taken into account through the introduction of void or vapor fraction which require 
empirical constitutive equations not always available In these much coarse grained approaches the 
detailed interface dynamics of bubbles and droplets is not available. 

There is a great recent interest in the computer simulations of interfacial dynamics for liquid- vapor co- 
existence. Prove of that is the large variety of techniques used to tackle the problem: Lattice-Boltzmann 
method Direct Montecarlo method Smoothed Particle Hydrodynamics (SPH) 10], Finite Differ- 
ence discretization of Navier-Stokes equations 0,0, and the Volume Of Fluid method Some of these 
techniques suffer from ad-hoc assumptions, more prominently isothermal behavior. Others require quite 
specific model technicalities in order to get sensible results or the specification of boundary conditions 
in complex topologies. It seems, therefore, appropriate to review here the theoretical foundation of the 
hydrodynamic equations for the van der Waals fluid. 

Although hydrodynamics equations for the van der Waals fluid were posed as early as 1901 by Kortcveg 

, | [ic| ] , a clarification of the structure and thermodynamic consistency of the equations is convenient in 
order to avoid potential problems. In the late 60's, for example, a variant of the theory was introduced 
by Kawasaki jl^, following the ideas of van Kampen [Q. It can be shown that Kawasaki theory is 
thermodynamically inconsistent. In Kawasaki theory, the long range attractive mean field potential 
between molecules acts as an external force in the momentum equation. A gradient approximation of 
this mean field potential leads to purely local equations involving third order spatial derivatives of the 
density field. Such a theory was proposed by Felderhof who constructed a thermodynamically consistent 
set of hydrodynamic equations for a van der Waals fiuid ||l^. However, Felderhof theory deals with the 
inviscid fiuid and no dissipation nor fiuctuations where included in his theory so a generalization of his 
theory seems appropriate. For an excellent review of diffuse interface models described by third order 
derivative terms see Ref. [ p^ . On the other hand, the discretization of third order derivatives is subject to 
high numerical errors. From a computational point of view, a treatment of the diffuse interfaces though 
a non-local integral term as in Kawasaki theory might represent an advantage worth exploring. 

The essential theoretical tool for describing equilibrium properties of liquid-vapor interfaces is the 
Density Functional Theory (DFT) ||l^ . A well-known local realization of the DFT known as the Density 
Gradient Theory successfully describes the equilibrium properties of liquid-vapor interfaces like surface 
tension and density profiles, and the correct relation between different scaling exponents near the critical 
point ijl^. This is due to the relatively smooth behavior of these interfaces as compared with fiuid-solid 
interfaces for which a non-local expression for the free energy functional is needed ||l^ . The work presented 
in this paper can be understood as a generalization of the local DFT to non- equilibrium situations in 
which an entropy functional instead of a free energy functional is used to describe both flow and thermal 
transport p7t|. 
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The strategy we follow is to derive the entropy functional from microscopic principles follo wing ideas 
of van Kampen, who derived a coarse grained free energy functional for the van der Waals gas flj]. This 
entropy functional is one of the building blocks of a recently developed framework for non-equilibrium 
thermodynamics known as GENERIC |18|] . This two generator formalism ensures thermodynamic consis- 
tency of its dynamic equations. Under a set of well-defined assumptions, the GENERIC formalism can be 
derived from first principles by means of a projection operator technique |p^ . It is therefore not surprising 
that all well-established dynamic equations for non-equilibrium systems fit into the GENERIC formal- 
ism: Relativistic and non-relativistic hydrodynamics, kinetic theory of gases (the Boltzmann equation) 
and of polymer systems in non-isothermal situations, chemical reactions, just to mention a few, have the 
GENERIC structure [Q. Significant new dynamical models for a number of complex systems, such as 
polymer/surface interactions ||T^, liquid crystals ||2^, or polymer blends ||2l[] have already been obtained, 
and even new cosmological models have been developed [ p2[ . In particular, the recognition of the role of 
the reversible dynamics in the system evolution equations has led to interesting new conclusions regarding 
closure approximations ||2^ and polymer reptation models |2j,|2^. It seems, therefore, sensible to formu- 
late a thermohydrodynamic theory for liquid-vapor coexistence according to the GENERIC framework. 
The basic bonus for this procedure is that strict respect for the first and second law of thermodynamics 
is guarded. 

The paper is distributed as follows. In Section || we present a brief summary of the GENERIC 
formalism and of the entropy functional for a van der Waals fluid. The entropy functional is computed 
explicitly in the Appendix [Vj. This appendix can be regarded as providing solid microscopic ground to 



the plausible assumptions made by Felderhof when postulating the entropy functional. In Section III 
we propose the hydrodynamic equations of a van der Waals fluid according to the GENERIC formalism 
and compare it with the theory of Kawasaki. In Section IV the local gradient theory is presented and 
compared with Felderhof 's theor y. So me c onclu ding remarks are given in Section while some further 



information is given in Appendix VIl , and VIII 



II. THE GENERIC FRAMEWORK AND THE ENTROPY 



A brief summary of the GENERIC structure is presented and we refer to the original references for 
further details [|l8| . The first essential step in the description of a system is the selection of the proper 
relevant variables x used to describe the system at a given level of description. The GENERIC dynamical 
equation for x is 

X = L{x)-VE{x) + M{x)-VS{x). (1) 

The first term in the right hand side produces the reversible part of the dynamics whereas the second 
term is responsible for the irreversible dissipative dynamics. Here, E{x), S{x) are the energy and entropy 
of the system expressed in terms of the variables x, V is the gradient operator in x-space, and L{x), M [x) 
are matrices that satisfy the following degeneracy requirements, 

L-V5 = 0, M-VE = Q. (2) 

In addition, L is antisymmetric (this guarantees that energy is conserved, i.e. the First Law) and satisfies 
the stringent Jacobi property ||l^. M is a positive definite symmetric matrix (this guarantees that the 
entropy is a nondecreasing function of time, i.e. the Second Law). If the system presents dynamical 
invariants I{x) different from the total energy, then further restrictions on the form of L and M are 
required 

\7I-L-VE:^0, V/-M-VS' = 0. (3) 

These conditions ensure that dl/dt = 0. The deterministic equations (|^) are, actually, an approximation 
in which thermal fluctuations are neglected. If thermal fluctuations are not neglected, the dynamics is 
described by a Fokker-Planck equation (FPE) that governs the probability distribution function p = 
p{x,t) [18|. This FPE has as equilibrium solution the Einstein's distribution function generalized to take 



into account the presence of dynamical invariants Eq, Iq |26|, this is 

p-^ix) = 6{E{x) - Eo)S{I{x) - lo) exp{S{x)/kB}, (4) 
where fc^ is Boltzmann's constant. 
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Our approach for constructing the entropy function for a hquid-vapor system is to compute micro- 
scopicahy p°'^{x) and, through Eqn. identify the entropy functionaL That is, we compute the joint 
probabiUty that an extended simple fluid in equihbrium has a particular realization of the mass, mo- 
mentum and internal energy density fields. We follow a cell-method first used by van Kampen in order 
to derive the equilibrium probability that a simple fluid has a particular realization of the mass density 



field The explicit details of the calculation are given in Appendix Vl. Physical space is divided 

in A4 cells of volume and the mass M^, momentum and internal energy e^j of cell fi in terms 
of microscopic coordinates of the molecules is considered. The internal energy contains the potential 
energy of interaction between particles within the same cell and the energy of interaction with particles 
of other cells. Provided that the microscopic coordinates are distributed according to the microcanonical 
ensemble we compute in Appendix [v| the joint probability P[A/, P,e] following similar steps as in Ref. 
||l2|. An essential step in the derivation is the assumption that the molecular potential has a repulsive 
hard core cj)^'^ (r) and a long range attractive tail 0' (r) . The interaction of particles between different cells 
is taken in mean field approximation and it involves only (fi{r). The final outcome is 



P[M, P, 6] = M, - 5 P,^ 



(5) 
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We recognize the dynamical invariants (total mass Mq, momentum Pg = and energy Eq) in the 
conserving delta functions. The entropy functional is given by 

S[M, e]^Y. ^'''(^A^ - ^M' ^'')- (6) 

Here, S^'^{e, M, V) is the entropy function of an isolated system oi N — M/m molecules (m is the mass 
of a molecule) interacting with only the hard core potential in a volume V. Note the non-trivial 

appearance of the mean field interaction energy cj) which is defined by 
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— Y,<P\R^,)M^M,, (7) 



where is the distance between cell centers and the prime denotes [i ^ v. This mean field potential 
involves only the long ranged part of the molecular potential. Integration of with respect to iV, e and 



using a steepest descent approximation (see Appendix VII) leads to van Kampen's expression, 

(8) 



P{M\ = exp {-/3F[M]} <5 - Mo^ , 



where the coarse-grained free energy functional is given by 

M 

F[M] = F'^'iP, M^,V^,) + 04>^. (9) 

Here, F^'^{f3,M^,Vfj) is the free energy of a system of hard core particles at temperature f3^^/kB = 
2K* / DM {K* is the most probable value of the kinetic energy in the system). Both Eqns. (^) and 
admit a continuum notation provided that the cells are large enough in order that the entropy and free 
energy are first order functions of their arguments 



^k,er] = J drs'^=(er-0r,Pr), (10) 
P[pr] - / dr(/-(pr,/3)+^r), (H) 



where er = e^/V^^, (j)^ = (pfj^/V^^ Pr = M^/Vfj,, and s^'^{er — <^riPr) is the entropy density of a hard core 
system, while f^^iPr, P) is the corresponding free energy density. The continuum form for the mean field 
potential energy (0) is 
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2m^ 



(12) 



The main benefit of the microscopic derivation presented in Appendix VII is the clear separation at 



the macroscopic level of the effects of the short and long ranged part of the microscopic potential, in 
particular the realization that the entropy of the full system is given directly in terms of the entropy of 
the hard core system, suitably evaluated at the corrected internal energy, Eqns. (||) or (|lO|). 

We now consider the particular situation in which the density field varies slowly in the scale of the 
range of 0'(r). In this case, we can Taylor expand pr' around r and retain up to second order in gradients 
in order to get 



-ap; 



(13) 



where we introduced a, c > through 



dr 



2to2 ' 



Substitution of ill'A) into (hlh leads to 



drr 



(14) 



F[Pr 



(15) 



which is the usual local gradient density expression for the free energy functional, also known as the 
Cahn-Hilliard or Ginzburg-Landau free energy functional for liquid surfaces psf . The corresponding 
entropy functional (p^, thus, has been computed at the same level of approximation. The coefficient c 
provides the overall magnitude of the surface tension coefficient [|5|. 
The energy function given by 



M 



2Mu 



(16) 



can also be written in continuum form as an energy functional 



E[pr,gr,er] = J 



2pr 



dr. 



(17) 



where gr = P^/V^ is the momentum density field. 



III. HYDRODYNAMIC EQUATIONS FOR A VAN DER WAALS FLUID 



In this section, we construct the hydrodynamic equations for a van der Waals fluid following the 
GENERIC formahsm. 

Due to the form in which the internal energy appears in the entropy, it proves convenient to use as the 
proper hydrodynamic variable the "intrinsic" internal energy el = — 4>j. which involves the interaction 
of molecules through the hard core part only. The energy (|7|) and the entropy (|l^) expressed in terms 
of these hydrodynamic variables are then 



E[p,g,e^] = 
S[p,g,e^] = 



1 

2 Pr 

s'^^ip.,el)dr. 



dr. 



(18) 



We have now the two basic building blocks for the GENERIC formulation of the dynamic equations 
(|l|) . The derivatives of the energy and entropy with respect to the relevant variables become functional 
derivatives with respect to the hydrodynamic variables, this is, 
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dE 

dx 



and 



dS_ 

dx 



SE 



ss 

<5gr 



2 r > 



(19) 



M''°(Pr,4) \ 



(20) 



The velocity field is Vr — gr/Pr, and the local temperature T^'^ and chemical potential /i^'^ per unit mass 
are 



he 



9s 



he 



de 



The 



ds 



he 



9pr 



(21) 



Note that these equations of state are those of the hard core system. By focusing on the reversible part of 
the dynamics, we must construct now the matrix L. A reasonable proposal is to use the same expression 
for L as in the case of hydrodynamics of a simple fluid. Note that if we neglect the interfacial potential 
energy (j)^ we should recover the hydrodynamics of a simple hard core fluid. We propose, therefore, for 
the matrix L — > Lrr' the following one 



/ Pr'V'Jrr' ^^ 

V ei,V'<5rr' + Pr'^'V'^rr' / 



(22) 



Here, V = d/dr' and 6rr' = S{r — r') is Dirac's delta function, P^'^ = P^^'^{pr, e*) is the pressure field as 
a function of the mass and intrinsic internal energy densities. 

As has been shown in Ref . , the matrix L satisfies the condition LVS = because of the Gibbs- 
Duhem relationship (for the hard core system). The fulfillment of LVS = is the basic reason for using 
P^'^ in Eqn. ( p^ ) instead of any other possible expression for the pressure. The matrix L also satisfies 
the stringent Jacobi identity With the explicit form ( |2^ ) for the L matrix and the gradient of the 
energy in Eqn. (19) the reversible part LVE of the equations of motion for the hydrodynamic variables 
are readily constructed 



dtPr 

dtSr 

dtei 



VVrgr - 



24 



Pr 

-Vej.Vr-P'^^(ej.,Pr)VVr. 



(23) 



Eqns. (g3|) are identical to the Euler equations of an inviscid (hard cored) fluid except for the additional 
term in the momentum equation that involves the long range part of the molecular potential. 

Concerning the irreversible part of the dynamics MVS, we observe that the same matrix M — > Mrr' 
that corresponds to the usual hydrodynamic equations is perfectly adequate. For simplicity, we assume 
that there is no extra dissipation due to the presence of interfaces. The matrix M is thus given by ||l8|| 



/ 
(W + lV-V')r?r7?'=5r 







Vo 



(24) 
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where 7 = Vvr + (Vvr)-^ is the symmetrized velocity gradient tensor. This matrix M satisfies the 
degeneracy MVE = as in the usual hydrodynamic case as the only modification of VE is in the 
term in the density component. We remark, however, that the transport coefficients (shear viscosity rjr-, 
bulk viscosity = kj. -\- 277r/3 and thermal conductivity Ar) that appear in M should be taken as state 
dependent and are expected to be strongly varying functions of the density field, in order to encompass 
the fact that vapor and liquid have very different transport properties. This is the reason for its space 
dependence displayed as a subindex r. The presence of T^^ in Eqn. ( p^ instead of other possible choices 
for the temperature, is dictated by the fact that MVS should produce the usual dissipative terms in the 
Navier-Stokes equation (note that VS in Eqn. ( pO| ) involves T^'^). 

The final set of hydrodynamic equations for the thermohydrodynamics of the van der Waals fiuid are 

dtPr = -VprVr, 

dtgr = -V-Vrgr - VP''%el,Pr) - V-T + p^F,, 

dtel = -Ve^Vr - P'^^el, Pr)V-Vr - V-q - r : Vvr, (25) 
where the constitutive equations for the viscous stress t and heat fiux are given by 

Tr - -Vr [VVr + Vv^] - (^r - 2r/r /3) ( V • Vr ) 1 , 

qr = -ArVTr'^^ (26) 



and the mean field force is given by 



-Vr^ /dr'0'(|r-r'|)pr' (27) 



Equations ( p5| ) are deterministic equations in which thermodynamic fluctuations are neglected. Ther- 
modynamic fluctuations are very easily included in the GENERIC formalism. One simply needs to take 
the square root in matrix sense of the M matrix, and this provides the amplitude of the noises , in 
accordance with the Fluctuation-Dissipation theorem. Because the matrix M in Eqn. ( |24| ) has the same 
structure of that of a one phase fluid, we can advance that the final fluctuating equations are obtained 
by simply adding to the stress tensor and heat flux qr a random counterparts given by p7[| 



fr = ^2T^^T]riar{t) - tr[ar(t)]l) + ^J 3T^^ nM^r (t)]! 

qr=T;^^V^Cr(i), (28) 

where fTr(i) is a matrix of delta correlated white noises and Cr{t) is a vector of delta correlated white 
noises (in space and time). Note that the form in which thermal fluctuations appear is similar to the 
approach proposed by Landau and Lifshitz However, one should note that in counterdistintion 

to their approach, here the random stress tensor and heat flux are multiplicative noises which are not 
Gaussian. What is Gaussian are the white noises (^r{t), Cr{t)- The strong dependence of the transport 
coefficients on the density due to the phase change has a direct effect on the amplitude of the thermal 
fluctuations. 

Associated to the stochastic differential equations describing hydrodynamic fluctuating variables there 



is a mathematical equivalent functional Fokker-Planck equation |27|. The GENERIC structure ensures 
that the equilibrium solution of this Fokker-Planck equation is given by the continuum version of Eqn. 
(^), as it should. 

Kawasaki proposed a set of hydrodynamic equations for a van dcr Waals fluid on intuitive grounds 
pT[ . Both the continuity equation and the momentum balance equation in ( p3| ) coincide with Kawasaki's 
postulated equations. On the other hand, he uses the entropy per unit mass — s^'^/pr instead of the 
intrinsic internal energy ej.. A standard calculation which uses the definitions (pi]) and the Euler equation 



yhcghc jj^^p — P^'^ — — shows that our equations (23) imply dts^'^ + Vr-Vs^'^ — 0. This is consistent 
with the fact that the reversible part of the dynamics should not produce an increase of the entropy in 
the system. In contrast, in Kawasaki treatment the reversible part of the substantial derivative of the 
entropy per unit mass is equated to prVr • V((/)r//Or) (in our notation). It is apparent that the addition 
of this term violates the conservation of the energy (llH) and should not be present for thermodynamic 
consistency. 
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IV. GRADIENT APPROXIMATION 



In this section we consider the gradient approximation ( p^ as apphed to the obtained equations (|25|). 
The last two terms of the momentum equation can be re-arranged in the form 

-VPi''' - PrV— = -VPr + CprVVVr, (29) 

where we have introduced 

Pr-P^^(Pr,4)-«P?- (30) 

The momentum equation thus becomes 

dtgr = -VVrgr - ^P{el, Pr) + VV^r- (31) 

It is iUuminating to use as the set of hydrodynamic variables Pr^^nU-r instead of Pr,gr,er- Here, the 
local internal energy is defined by 

Wr = ej. — ap^, (32) 

which, from a microscopic point of view, represents the potential energy due to local interactions, that 
is, hard core interactions plus the mean field attraction "within" the small volume at r. By taking the 
time derivative of ( ^2|) and using the hydrodynamic equations (|2^), a standard calculation leads to 

dtPv = -VprVr, 

a^gr - -VVrgr - VP"'^W(Ur, Pr) + CprVV^Pr - V • T, 

^^U,, = -VWrVr - P"'^^(Ur, Pr)V-Vr - V-q - T : VVr, (33) 

where the heat fiux is given by 

qr = ~ArVr^<^W(Mr,Pr). (34) 



Here, pvdWjvdw ^-^^ equations of state for the va n de r Waals fiuid (see App endix |VIII[ ). In 



obtaining (^) and (p4[), use has been made of the identities (IOC) in the Appendix |VII]|. 

The hydrodynamic equations (^3|) are equivalent to those proposed by Felderhof |13| . However, in 
Felderhof 's theory, dissipation was neglected and, consistently, thermal fiuctuations where not considered. 
Thus, equations ( P3|) with (|2^) and (Q) can be taken as the natural generalization of Felderhof theory. 

It is not obvious at a first glance that the total momentum of the system is conserved by the term 
cprVV^Pr in (^3|). That this is the case can be seen by noting that this term can be cast into the 
form — V-Pr, and the momentum equation has the form of a balance equation. The "surface tension" 
contribution to the stress tensor Pv has the form EOtl 



2c 



ivprVpr + -^(Vpr)'l ~ ^PrVVpr ^ ^PrV^Prl 

6 12 3 6 



(35) 



This contribution to the stress tensor is identical to the one derived from molecular considerations for 
the equilibrium stress tensor [ pO[ under the gradient approximation p^ . Wc observe, therefore, that the 
theory that is being presented relics also on the hypothesis of local equilibrium, in much the same way as 
in the usual continuum hydrodynamics approach. 



V. SUMMARY AND DISCUSSION 



We have computed the entropy functional for a simple fluid starting from microscopic principles under 
the classic van der Waals assumption that the molecular potential has two well-separated ranges, a 
hard core repulsive part plus a very long ranged attractive tail. Under the further assumption that the 
hydrodynamic fields vary slowly on the total range of the potential it is possible to write a local entropy 
functional that depends on the density gradients. We have also presented the connection between this 
entropy functional and the local gradient free energy functional which is used in the context of the 
study of equilibrium of liquid-vapor interfaces. Even though the assumption of the unrealistic molecular 
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potential seems to be too restrictive, the equilibrium gradient theory is highly successful in describing 
the equilibrium properties of these interfaces, and gives confidence on the dynamical approach taken in 
this paper. 

The entropy functional describes the equilibrium properties of the system but it is also one of the basic 
inputs for constructing the non-equilibrium evolution equations for the hydrodynamic fields through the 
GENERIC formalism. By making use of all the well-known information about the L and M matrices for 
the hydrodynamics of a single phase fluid flow, we can construct the hydrodynamic equations for a van 
der Waals fluid in a rather simple way. 

The hydrodynamic equations proposed predict unstable forces if the thermodynamic state is in the 
positive slope part of the van der Waals pressure diagram. When this happens, the fluid spontaneously 
adopts one of the two possible densities corresponding to liquid and vapor. Density variations appear in 
the short length scales associated to the interfaces of bubbles or drops. It is precisely in these interfacial 
regions where the last term in Eqn. is important. This term generates forces in the fluid associated to 
the surface tension of the interface. Actually, the interfaces try to adopt a spherical shape. The interface 
is a diffuse object with a finite width. Note that no boundary conditions are required on the interface of 
a bubble or a drop. Actually, there is no need to known a priori the location of the "interface" in order 
to have a well-posed hydrodynamic problem. 

The set of equations (|2^) with the gradient approximation (^3|) and the set of equations ( ^3| ) are 
mathematically equivalent. We regard as the main benefit of the microscopic derivation presented in this 
paper the clear physical interpretation of the variables (either or Mr) and the correct equations of state 
(either P^"^ or p^^w^ j-,g ^ggjj jjj ^j^g equations. This issue is not always clear in the usual presentations 
of this subject based on phenomenological arguments that simply add gradient terms either in the energy 



From the computational point of view of numerically solving the hydrodynamic equations obtained, we 
observe that a theory like Kawasaki's in which the mean field long range forces appear in integral form, 
Eqns. (^5|), might represent an advantage in front of a theory like Felderhof's in which these forces appear 
in differential form through third derivatives of the density field, Eqns. (p3|). In those regions where the 
interfaces are located, the density field changes strongly which means that the third derivatives present a 
rather complex structure in very short length scales. In order to get stable and convergent results, very 
fine grids are required. These problems do not arise in the discretization of the integral form of the long 
range forces. 

The theory derived here connects with the well-known phase field theories that deal with systems 
that present complex boundary conditions like those occurring in melting and dendritic growth pl| | or 
immiscible fluids ||3^ . In these systems, the partial differential equations governing the evolution of the 
relevant variables require the formulation of boundary conditions on moving surfaces. The motion of the 
surface is, in turn, coupled in a complex way with the dynamics of the relevant variables. The approach 
taken by the phase field method is to include the boundary conditions into the equations of motion 
through the introduction of an auxiliary phase field that locates in a somewhat diffuse way the position 
of the interface. A further dynamics is generated for the phase field which is coupled with the dynamics 
of the rest of fields in the system. This dynamics is obtained from a phenomenological formulation of a 
"free energy". The theory we have derived is then in the class of a phase field theory for bubbles and 
drops. In our case, the phase field has a definite physical meaning in terms of the density field, and no 
phenomenological free-energy needs to be introduced ad hoc. 

P.E. would like to thank useful discussions with H.C. Ottinger, P. Tarazona, J. Cuesta, and R. Delgado. 
This work has been partially supported by DGYCIT PB97-0077. 

VI. APPENDIX: MICROSCOPIC DERIVATION OF THE ENTROPY FUNCTIONAL 

In this section, we compute from first principles the joint probability that an extended simple fluid in 
equilibrium has a particular realization of the mass, momentum and internal energy density fields. This 
generalizes van Kampen derivation which only considered the mass density field |12|] . 

The isolated fiuid is assumed to be described microscopically by 2;, the set of positions q; and momenta 
Pi of the center of mass of its Aq constituent molecules (no rotational or vibrational degrees of freedom 
are considered, for simplicity). The Hamiltonian of the system is given by 



or the entropy [|4j. 




(36) 
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where m is the mass of a molecule and 0y is the pair- wise potential function between molecules i and j, 
assumed to depend only on its separation . The total momentum of the system is given by 

No 
i 

The microscopic dynamics is such that H{z) and P(2) are dynamical invariants. Obviously, the total 
number A'o of particles (or the total mass) is also a dynamical invariant. 

On a mesoscopic scale, the space containing the fluid is partitioned in M cells. The state of the system 
at that scale is described with the mesoscopic variables N^, P^, e^, /i = 1, . . . , M, where is the number 
of particles in cell /i, is the momentum of the center of mass of the particles of cell fi and e^i is the 
energy with respect to the center of mass of the particles in cell ^. In order to relate the mesoscopic 
variables with the microscopic ones, we introduce the characteristic function XmI'") '^^ M which takes 
the value 1 if r is within cell /i and otherwise. For example, the space can be partitioned into Voronoi 
cells whose characteristic function is Xm('^) = H^^^ ^(k ^ '^A ^ |r — R/^l) where 9{x) is the Heaviside 
step function and R^^ is the center of the Voronoi cell /i. Obviously, 

M 

because the cells cover all space without overlapping. The volume of the cell is given by 

V^.= I d'rx^(r), (39) 

JVa 

where Vb is the total volume of the system. The mesoscopic variables N^^ P^ and can now be written 
as functions of the microstate z, this is 



i 
i 

vw-E(i(p.-^)^*.) 



xMi)^ (40) 
where we have introduced the potential energy of particle i through 
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(41) 



Note that we use Latin indices to refer to variables at the microscopic level and Greek indices to refer to 
mesoscopic variables. 

The dynamical invariants can be written in terms of the mesoscopic variables by using Eqn. ( ^Sf ) , this 

is 

^(^) = E^ + ^.(^)^ 

V{z)^Y.^,{zl (42) 

where AI^ = mN^ is the mass of cell p,. 

We will assume that the system is at equilibrium, which means that the effective probability density 
p°'^{z) of finding a particular value of z is a function of the dynamical invariants of the system only. If 
we assume that these invariants are known with precision, the equilibrium distribution function is given 
by the molecular ensemble 

P'^'i^) = ]^'^(P(^) - ^^W{z) ~ Eo), (43) 



9 



where Pq is the total momentum of the system and Eq the total energy. The factor Nq] is quantum 
mechanically in origin (it comes from the undistinguishability of the particles) and solves the Gibbs 
paradox (and makes the macroscopic entropy a first order function of energy and number of particles). 
The normalization factor Qq is given by 

f dz 

Qo = niUo, No, Vo) = J j;^SiP{z) - Po)<5(i/(^) - Eo). (44) 

Here, Uq — Eq — Pq/2Mq is the internal energy of the total system and the total volume Vb appears 
parametrically as the region of integration of positions. The macroscopic entropy for this system is defined 
as 

^mac(^^ TV, V) = kB In r!o([/, V). (45) 
The intensive parameters are defined as the derivatives of the macroscopic entropy, this is, 



1 




T 


dU 


M 




T 


dN 


P 


^^mac 


T 


dV 



(46) 



where T is the temperature, ^ the chemical potential, and P the pressure. 



A. Calculation of P[N, P, e] 

The probability that a set of functions X(z) take particular values x when the system is at equilibrium 
is given by 

(47) 

where p'^'^{z) is the equilibrium ensemble. We now consider the functions X{z) that appear in Eqn. ( ^7|) 
to be the set of mesoscopic variables (pO|). In this way, the probability that the system adopts a particular 
set of values iV^, P^, for each cell is given by 

M 



P[iV,P,e]= j dzp-'^{z)f{x{N^{z)~N^) 



X 5(P^(z)-P^)5(e^(z)-e^). (48) 

In this expression, we have introduced the characteristic function x {Nf_i{z) — iV^) that takes the value 1 
(not infinity) for the region of phase space in which the microstate z produces exactly the value for 
the number of particles in cell /x. The probability (^8|) is normalized according to 

„ „ Wo 

dPi...dPM dei...deM J2 PiN,'P,^] = l: (49) 

Ni,...Mm 

where the numbers Nfj, are subject to J^fi-^i^ ~ -^o- By using the equilibrium ensemble ( |4^ ) and Eqns. 
(E^) we obtain 



/ M 



/ ^J[xiN,iz)-N,)SiP,{z)-P,) 
X ,5(e^(z)-e^). (50) 



10 



The term Yl^ x{^i-i{^) ~ within the integral is zero unless the microstate z is such that there are 
exactly Ni particles in the first cell, N2 in the second, etc. There are NqI/NiI ■ ■ ■ Nm^- ways of having the 
No particles distributed among the M cells with the prescribed numbers in each cell. Therefore, we can 
write 

M 




M / N,, 



J dpi... dpM f[S Ph' ~ ^ 




(51) 



where the function ^{P, E, N) is introduced and computed in the Appendix VII. The numbers 
Ni, . . . , Nm are subject to J2fj,^f^ ~ -^0 and, therefore, we can extend the above result to arbitrary 
iVi, . . . , Nm provided that we multiply Eqn. ( Jsi] ) with the factor xij^^ ^fj. ~ ^a), which takes the value 
1 if its argument is zero. 

Let us write now the potential energy of cell /z as follows, 

E^v = ^EE^w. + ^EEE^«.- (52) 

The first term is the potential energy due to the particles of cell /i and it is a function solely of the 
coordinates of the particles that are in cell jj,. The second term, the potential energy due to the interaction 
of particles in different cells, involve coordinates of particles in cells v different from ^i. It is this last term 
which hinders the decoupling of the integrals in positions in Eqn. (pT]). For this reason, it is necessary 
to make approximations to this last term. In what follows, we will approximate the potential energy due 
to different cells by a mean field constant (independent of coordinates of particles), this is, 

^ E E E -^v.. « E ^^M^^c-M- (53) 

This approximation is fully justified when the potential (\)(r) has two typical ranges, a short length scale 
Ts much smaller than the typical size of Voronoi cells and a long length scale r; much larger than the size 
of the cells. In this case, we can separate 0(r) = 0*(r) + (/''(r) and write 

^EE^«.-^EE<..' (54) 

because the short range part is negligible for different cells ji ^ v. Also, because the range r; is much 
larger than the size of a cell, we can further approximate Eqn. (^) as if all particles in each cell where 
located at the center of the cell and, therefore, we can write 

2 E E ^ -4>\R^.)N^.N,, (55) 

V 3" 
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and the potential energy due to different cells has the structure of Eqn. ( pc_ ^ 

Under the mean field approximation (|5^), the different position integrals in ( ^l|) decouple, so that Eqn. 
(pWj can be written as 



M 

n 



(56) 



where 



(57) 



where i/)^ is the mean field inter-cell potential energy of cell /j. 
We can write Eqn. ( p6|) as 

M 



(58) 



where ^^^'^{U, N, V) is defined in Eqn. (44). Note that the Hamiltonian that should appear in Eqn. ( [44|) 
is that of a system of particles interacting through the short range part of the potential (t)^{r). Finally, 
Eqn. (BG) can be written as 



M 



ME 



2Af„ 



^0 



X exp{5[7V,e]/fcB}, 



where the entropy functional at the mesoscopic level is defined by 

5[A^,e]=^5'-(e,,-0^,iV,„V^), 



(59) 



(60) 



where we have used Eqn. (|45|). The mesoscopic entropy is given in terms of the sum of the macroscopic 
entropies of each cell, as if they were isolated systems in a volume V^, with number of particles N^j_ and 
with an energy e^ — (f>^. The internal energy in (^) includes in its microscopic definition the potential 
energy of particles that interact with particles of neighboring cells, that is, includes the inter-cell 



energy. In this way, 



represents the purely internal energy due to the particles within the cell (in 



mean field). Even though we have made progress in writing the probability for the mesoscopic variables, 
Eqn. (p9|), the entropy of each cell S''^'^(e^, V^, iV^) is still an unknown quantity. We remark that this 



entropy is the macroscopic entropy of a system of molecules interacting with the short range part of the 
potential. In applications of the model, this entropy will be given by simple models (see Appendix VIII ) 
or by very accurate expressions like the Carnahan-Starling equation of state p^ . 

Apparently, Eqn. ( |60|) simply says that "the entropy is additive" . However, this sentence is imprecise: 
The entropy S[N, e] in the Ihs of Eqn. (|6^) is a different object from the entropy S'^'^{e^, N^, V^) appearing 
in the rhs. They depend on a different number of variables. Therefore, one should rather say that "the 
mesoscopic entropy at the level of hydrodynamic variables is the sum of the entropy at the level of 
dynamical invariants (the macroscopic entropy) of each cell as if they were isolated" . 

The macroscopic entropy is a first order function of its variables and therefore. 



S[N, 



(61) 
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which admits the continuum notation 

5K,er] = / drshc(er_0^,nr), (62) 



where we have introduced the continuum version of the effective interaction energy between molecules in 
different points of space (cells), 

T.^,^ I dr, 
V ' 



(63) 



Note that in this continuum notation we can write the dynamical invariants ( [4 21 ) of the system as 

fdrfl- 



^ 2 mrir 

P - / drgr, (64) 



where P^/V^ gr is the momentum density field. We stress that we regard Eqn. (62) not as an strict 
mathematical limit in which the volume of the cells tend to zero, but rather as a convenient notational 
tool. 

In spite of its apparent form, the mesoscopic entropy ( |6^ ) is not local in space. This is, it is not given 
by a sum of a function of the variables at a given point. This is due to the presence of the term (p^ which 
can be written in continuous notation as 



(j)^= J dr'(j)'{\r-r'\)nrnr'. (65) 

We have implicitly associated the potential energy of interaction with other cells (p^. with the phenomena 
of surface tension. In order to make this connection more explicit, in the following subsections we consider 
the marginal probability distribution functions and will make contact with the Density Functional Theory. 



B. The marginal distribution P[N,e] 



By integrating the distribution function P[N,F,e\ over momenta we will have the probability of a 
realization of the "fields" N, e irrespective of the values of the momenta in each cell, this is 



P[N, e] = exp {S[N, eykg} - N^^ 

X JdP,...dPMS (Y.P^^-Po^ 



M 



2M^ ^ 



-En 



exp {S[N, e]/kB} ^X - No^ 



M 



D(M-l) 



(66) 
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where we have used once more Eqn. (^0|) in the appendix. 

The most probable realization iVj*, . . . , N^j, el, ... , e*j^j of the fields is the one which maximizes the 
functional 



D 



(67) 



where we have introduce the Lagrange multiplier (3X that takes into account the restriction A^^ — Nq 
The maximum N* , e* of the functional ( |6^ ) is the solution of the following set of equations 

dS 



dS 
de,, 



[N*,e*] =kBf3X, 



D{M-1) 



(68) 



We have neglected a term of order N^^ in the first equation, which is reasonable if typically there are 
many particles in each cell. Eqns. ( pq ) is a set of 2M + 1 equations for the unknowns N*,e*^,X. Note 
that the solution TV*, e* , f3X depends parametrically oiq A'o, Uq. 



We find now a convenient approximation to Eqn. (|66|) by noting that this probability is expected to be 
highly peaked around the most probable state. Therefore, for those values of the field e for which P[N, e] 
is appreciably different from zero (that is, around e*), we can approximate 



P 



exp{/? E 



0} 



ctn.exp{-/3E^A'}' 



(69) 



where P = D{M — l)/2 — 1 is a very large number and wc have introduced 

D{M -l)/2-l 



(70) 



Note that /3 is proportional to the inverse of the most probable value of the kinetic energy per cell. 
Finally, we can write Eqn. (|66|) as 



P[N, e] = ^ exp 1 5[7V, e] /fcs - /3 E I 



(71) 



where is the corresponding normalization function. We see, therefore, that by integrating the momenta 
the "microcanonical" form Eqn. (Mh becomes the "canonical" form (O). 



14 



C. The marginal distribution P[N] 



There are two different routes to compute the probabihty of a certain distribution iVi , . . . , N]\j of the 
particles in the cehs. We could simply start the calculation from Eqn. ( ^8| ) without the momentum and 
energy conserving delta functions. This route is essentially the one taken by van Kampen [Q. The 
alternative is to integrate out the energy field in Eqn. @. The same results are obtained in both 
approaches and we illustrate here this second one, this is, 



P[N] - J dei... deMP[N, e] 

= exp |-/3^ Fh^(/3, N^, V^) - /3^^| 

^ M \ 



^X\^2^N^-N„j, (72) 

where we have used the form ( pO| ) for the entropy function, we have changed variables to = — 4>^ 
and have introduced the macroscopic free energy F^'^{/3, N,V) of a system of N particles at inverse 
temperature /3 in a volume V, and interacting with a short range potential. In deriving Eqn. (^2|) we 
have made use of the following identity 

poo 

/ eiip{S{U)/kB - l3U}dU = exp{-PF{P)}. (73) 
Jo 

This expression is the statistical mechanics link between the macroscopic entropy and the free energy 
and can be proved by computing the Laplace transform of fla in Eqn. (H). One has 

oo 

no{U,N, V)exp-{l3U}dU 

= f dzS{J2p,)exp-{pHiz)} = Zip.N -1,V), (74) 
i 

where the partition function Z{P) is defined through this equation. The presence or absence of the 
momentum conservation delta function is irrelevant when the number of particles is large so we may very 
well drop it. The free energy F(/3) is defined as 

/3F(/3) = -lnZ(/3), (75) 

The usual thermodynamic link between entropy and free energy, 

F{T,N,V) = U -TS{U,N,V), (76) 

can be obtained under the assumption that the integrand in ([73| ) is highly peaked. 

VII. APPENDIX: MOLECULAR ENSEMBLE 

In this appendix we compute explicitly the following integral 

ci>(Po, i?o, M) = J d^^'p ^ f E ^ - ^0 j 

x<5^(^EP'-Po)' (77) 
which appears repeatedly when computing molecular averages. By a simple change of variables we obtain 
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M / M \ 

$(Po, £^0, M) ^ n(2m,)^/2 j d^^^p 5 ( ^ p2 - j 

X 5^|^^(2m.)i/2p,_Poj , (78) 
The equation ^^^(2m.i)^/^pj = Pq are actually D equations (one for each component of the momentum) 



which define D planes in R . The integral in (|77|) is actually over a submanifold which is the intersection 

of the D planes with the surface of a DM dimensional sphere of radius E^'^ . This intersection will be 
also a sphere, which will be now of smaller radius and also of smaller dimension, D{M — 1). 
In order to compute (|77|), we change to the following notation 

V ^ (pT, • • ■ ,PM,Fi, ■ ■ ■ ,Pm,p1, ■ ■ ■ ,Pm) 

= ((2mi)i/2,...,(2mM)'/',0,...,0,0,...,0) 
Cy^{0,...,0, (2mi)i/2, . . . , (2mM)'/', 0, . . . , 0) 

= (0, . . . , 0, 0, . . . , 0, (2mi)i/2, . . . , (2mM)'/'). (79) 

Note that these vectors satisfy Cx Cy — Cy Cz — Cz Cx = 0. With these vectors so defined, Eqn. 
( |77| ) becomes 



X 5{CyV-P^) 5{Cz-V~P§) 



'V 

M 



X 



\{{2m,)'"^. (80) 



Now we consider the following change of variables 

C C 



r' = v- -^ps + jp^p^ + jp^Po' , (81) 



\Cx\' ° " |C,|2 ^ |C.|2 

which is simply a translation and has unit Jacobian. Simple algebra leads to 

M . 

<i>(Po,£;o,M) = Y[{2m.f/' / d^^'r' 6{V'^ - [/q) 
i 

X 5{Cx-V') S{CyV') S{Cz-r'), (82) 
where we have introduced the total internal energy 

where Mo = J2i "^j *^he total mass. 

We now consider a second change of variables V" = A-T'' through a rotation A such that 



ACx 


= (1,., 


.,0,0,. 


,.,0,0,. 


..,0) 


ACy 




.,0,1,. 


,.,0,0,. 


..,0) 


AC, 


= (0,.. 


.,0,0,. 


,.,0,1,. 


..,0) 



(84) 

It is always possible to find a matrix A that satisfies Eqns. (^J). For example, consider a block diagonal 
matrix made of three identical blocks of size M x AI. Then assume that each block is the same orthogonal 
matrix which transforms the vector ((2mi)^/^, (2m2)^/^, . . . , (2mM)^^^) into (2A^o)"^''^(l, 0, . . . , 0). After 
the rotation (which has unit Jacobian and leaves the modulus of a vector invariant) the integral ( ^^ 
becomes 
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M 

$ = [](2m,)^/2 / d''^'V"S{r"^ - Uo) 



X sip'n 6{p';y) s{pr) 

M 



i 

We compute now the integral over the sphere in Eqn. ( |85| ) by using that the integral of an arbitrary 
function F(x) = /(|x|) that depends on x only through its modulus |x| can be computed by changing to 
polar coordinates 

j F{^)d^'^^ujM f{r)r^^-^dr. (86) 

The numerical factor ljm, which comes from the integration of the angles, can be computed by considering 
the special case when /(r) is a Gaussian. The result is 

^M/2 



By using Eqn. (86), Eqn. (|85| ) becomes 

M 

i 

X j dp pD(M-l)-l S^p2 _ 

We need now the property 

where Xi are the zeros of f{xi) = 0. For the case of Eqn. ( p8| ) we have f{x) = x"^ — U , Xi = ±(C/)^/^ and 
f'{x) = 2x. Therefore, 

1 D(M~1)~2 

<!>{Po,Eo,M)^ -u;DiM-i)Uo ' [|(2m,)^/2. (90) 



VIII. APPENDIX: VAN DER WAALS AND HARD CORE MODELS 

The particular functional form of s^^{e, p) cannot be computed from the microscopic analysis presented 
in the previous appendix. The approach taken in this paper is that this fundamental thermodynamic 
equation (in the sense of Callen |]33[|) is known either from empirical sources or by suitable modelling. In 
this appendix we summarize the results for the fundamental equation of van der Waals model and for its 
corresponding hard core model. 

The van der Waals model can be defined through the following fundamental equation, which relates 
the entropy density with the internal energy density e and number density n, 

s "(e, n) = — - — fcsn-fcsnlnl-t ^ _ ]■ (91) 

The number density is n = p/rriQ where p is the mass density and toq is the mass of a molecule. The 
thermal wavelength is defined by 

AvdW(^^ „) = (2^^^^^^vdW(g^„))i/2' (92) 
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and we have introduced the function 

The two equations of state are the derivatives of the entropy with respect to each variable 



(94) 



where r^^* is the temperature and the chemical potential. Straightforward calculations lead to 

.-(.n) = *,r-" (,„ + ^) - .an. (95) 

The third equation of state can be obtained from the Euler equation P = Ts — e + /xn, which leads to 

P (e,n)=n ^ — -—an. (96) 

1 — no 

The hard core model that corresponds to the van der Waals model can be defined as the model that 
results from taking the attractive parameter a = in the van der Waals model. This model is essentially 
the gas ideal model but with excluded volume effects. In this way, one obtains 

s {e, n) = -^—ksn - ken In ^ ^ \' , (97) 



1-nb 



where now 



A^"^(e,n) = 



(27rTOoA:BTh<=(e,n))V2^ 



T^^(e,n) = -^-^. (98) 
D ksn 



The equations of state are now 

T'^^(e,n)=T^^(e,n), 



1 — nb J 1 — nb 

P-^ie,n)=n'-^^^'. (99) 

Note that the following functional relations hold 

T^^u + an^,n) =T'''^'^{u,n), 
P^'iu + an^, n) - an^ = P"'*^(m, n), 

l?%u + an'^, n) - 2an = n^^'^iu, n). (100) 
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